/*
 * GaseousActivity.cpp
 *
 *  Created on: 20 Aug 2011
 *      Author: Allan
 */

#include "GaseousActivity.h"

// GeoReact includes
#include <Geochemistry/Activity/GaseousModels/SpycherReedCO2Model.h>
#include <Utils/Algorithms.h>

GaseousActivity::GaseousActivity() 
{}

GaseousActivity::GaseousActivity(const string& species, const vector<string>& phase_species, const ActivityCoefficientFunction& activity_coefficient) :
activity_coefficient(activity_coefficient)
{
	index_species = IndexOf(species, phase_species);
}

GaseousActivity::~GaseousActivity()
{}

const double GaseousActivity::activity(double T, double P, const VectorXd& n) const
{
	const double molar_fraction = n[index_species]/n.sum();
	
	return P * molar_fraction * activity_coefficient(T, P, n);
}

const double GaseousActivity::activityCoefficient(double T, double P, const VectorXd& n) const
{
	return activity_coefficient(T, P, n);
}

const VectorXd GaseousActivity::activityDerivatives(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	const unsigned num_gaseous_species = n.rows();
	
	const double n_total = n.sum();
	
	const double average = n_total/num_gaseous_species;
	
	// Check if all gaseous species have the same number of moles (this is to avoid singular matrix)
	const bool same_number_moles = ((n.array() - average)/average).abs().sum() < 1.0E-6;
	
	const double auxiliary_value = (same_number_moles) ? 0.0 : -a[index_species]/n_total;
	
	VectorXd derivatives = VectorXd::Constant(num_gaseous_species, auxiliary_value);
	
	derivatives[index_species] += a[index_species]/n[index_species];
	
	return derivatives;
}

const GaseousActivity GaseousActivityIdealModel(const string& species, const vector<string>& phase_species)
{
	// The ideal model for the gaseous activity coefficient (gamma = 1.0)
	auto ideal_gaseous_activity_coefficient = [](double T, double P, const VectorXd& n) { return 1.0; };
	
	GaseousActivity gaseous_activity(species, phase_species, ideal_gaseous_activity_coefficient);
	
	return gaseous_activity;
}

const GaseousActivity GaseousActivitySpycherReedCO2Model(const vector<string>& phase_species)
{
	GaseousActivity gaseous_activity("CO2(g)", phase_species, SpycherReedFugacityCoefficientCO2);
	
	return gaseous_activity;
}

const GaseousActivity GaseousActivitySpycherReedH2OModel(const vector<string>& phase_species)
{
	GaseousActivity gaseous_activity("H2O(g)", phase_species, SpycherReedFugacityCoefficientH2O);
	
	return gaseous_activity;
}
